﻿using System;
using System.Text;
using SoftwareConsulting.BI11.NumericalMethods.Iterations;

namespace SoftwareConsulting.BI11.NumericalMethods.MatrixAlgebra
{
    /// <summary>
    /// Класс метода Якоби для нахождения собственных значений и векторов матрицы
    /// </summary>
    public class JacobiЕigenValuesAndVectorsFinder : IterativeProcess
    {
        private double[,] _rows;
        private double[,] _transform;
        private int _p, _q;	//индексы самого большого недиагонального элемента

        /// <summary>
        /// Конструктор класса
        /// </summary>
        /// <param name="simmetricMatrix"></param>
        public JacobiЕigenValuesAndVectorsFinder(SymmetricMatrix simmetricMatrix)
        {
            int n = simmetricMatrix.RowsCount;
            _rows = new double[n, n];
            for (int i = 0; i < n; i++)
            {
                for (int j = 0; j < n; j++)
                    _rows[i, j] = simmetricMatrix[i, j];
            }
        }

        /// <summary>
        /// Возвращает собственные значения матрицы
        /// </summary>
        /// <returns></returns>
        public double[] GetEigenValues()
        {
            int n = _rows.GetLength(0);
            double[] eigenvalues = new double[n];
            for (int i = 0; i < n; i++)
                eigenvalues[i] = _rows[i, i];
            return eigenvalues;
        }

        /// <summary>
        /// Возвращает собственные вектора матрицы
        /// </summary>
        /// <returns></returns>
        public Vector[] GetEigenVectors()
        {
            int n = _rows.GetLength(0);
            Vector[] result = new Vector[n];
            double[] temp = new double[n];
            for (int i = 0; i < n; i++)
            {
                for (int j = 0; j < n; j++)
                    temp[j] = _transform[j, i];
                result[i] = new Vector(temp);
            }
            return result;
        }

        /// <summary>
        /// Выполняет итерацию
        /// </summary>
        /// <returns></returns>
        public override double DoIteration()
        {            
            double nonDiagonal = GetLargestNonDiagonalValue();
            Transform();
            return nonDiagonal;
        }

        /// <summary>
        /// Выполняет обмен элементов согласно методу Якоби
        /// </summary>
        /// <param name="rowIndex"></param>
        private void Exchange(int rowIndex)
        {
            int row2Index = rowIndex + 1;
            double temp = _rows[rowIndex, rowIndex];
            _rows[rowIndex, rowIndex] = _rows[row2Index, row2Index];
            _rows[row2Index, row2Index] = temp;
            int n = _rows.GetLength(0);
            for (int i = 0; i < n; i++)
            {
                temp = _transform[i, rowIndex];
                _transform[i, rowIndex] = _transform[i, row2Index];
                _transform[i, row2Index] = temp;
            }
        }

        /// <summary>
        /// Завершает итерации
        /// </summary>
        public override void FinalizeIterations()
        {            
            int n = _rows.GetLength(0);
            int bound = n - 1;            
            while (bound >= 0)
            {
                int m = -1;
                for (int i = 0; i < bound; i++)
                {
                    if (Math.Abs(_rows[i, i]) < Math.Abs(_rows[i + 1, i + 1]))
                    {
                        Exchange(i);
                        m = i;
                    }
                }
                bound = m;
            }
        }

        /// <summary>
        /// Подготавливает итерации
        /// </summary>
        public override void CreateRandomPopulation()
        {
            _transform = SymmetricMatrix.IdentityMatrix(_rows.GetLength(0)).Components;
        }

        /// <summary>
        /// Возвращает самый большой недиагональный элемент
        /// </summary>
        /// <returns></returns>
        private double GetLargestNonDiagonalValue()
        {
            double result = 0;            
            int n = _rows.GetLength(0);
            double r = 0;
            for (int i = 0; i < n; i++)
            {
                for (int j = 0; j < i; j++)
                {
                    r = Math.Abs(_rows[i, j]);
                    if (r > result)
                    {
                        result = r;
                        _p = i;
                        _q = j;
                    }
                }
            }
            return result;
        }

        /// <summary>
        /// Возвращает матрицу как строку
        /// </summary>
        /// <returns></returns>
        public override string ToString()
        {
            StringBuilder sb = new StringBuilder();
            char[] separator = { '[', ' ' };
            int n = _rows.GetLength(0);
            for (int i = 0; i < n; i++)
            {
                separator[0] = '{';
                for (int j = 0; j <= i; j++)
                {
                    sb.Append(separator);
                    sb.Append(_rows[i, j]);
                    if (j == 0) 
                        separator[0] = ' ';
                }
                sb.Append('}');
                sb.Append('\n');
            }
            return sb.ToString();
        }

        /// <summary>
        /// Преобразует матрицу по методу Якоби
        /// </summary>
        private void Transform()
        {
            double apq = _rows[_p, _q];
            if (apq == 0)
                return;
            double app = _rows[_p, _p];
            double aqq = _rows[_q, _q];
            double arp = (aqq - app) * 0.5 / apq;
            double t = arp > 0 ? 1 / (Math.Sqrt(arp * arp + 1) + arp)
                               : 1 / (arp - Math.Sqrt(arp * arp + 1));
            double c = 1 / Math.Sqrt(t * t + 1);
            double s = t * c;
            double tau = s / (1 + c);
            _rows[_p, _p] = app - t * apq;
            _rows[_q, _q] = aqq + t * apq;
            _rows[_p, _q] = 0;
            _rows[_q, _p] = 0;
            int n = _rows.GetLength(0);
            for (int i = 0; i < n; i++)
            {
                if (i != _p && i != _q)
                {
                    _rows[_p, i] = _rows[i, _p] - s * (_rows[i, _q] + tau * _rows[i, _p]);
                    _rows[_q, i] = _rows[i, _q] + s * (_rows[i, _p] - tau * _rows[i, _q]);
                    _rows[i, _p] = _rows[_p, i];
                    _rows[i, _q] = _rows[_q, i];
                }
                arp = _transform[i, _p];
                aqq = _transform[i, _q];
                _transform[i, _p] = arp - s * (aqq + tau * arp);
                _transform[i, _q] = aqq + s * (arp - tau * aqq);
            }
        }
    }
}
